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Abstract: We discuss the optimal detection of point sources from multiwavelength imaging data using 
an approach, referred to as MDET, which requires no prior knowledge of the source spectrum. MDET 
may be regarded as a somewhat more general version of the so-called "chi squared" technique. We 
describe the theoretical basis of the technique, and show examples of its performance with four-channel 
infrared broad-band imaging data from the WISE mission. We also discuss the potential benefits of 
applying it to the multifrequency data cubes of the ASKAP surveys, and suggest that it could increase 
the detection sensitivity of searches for neutral hydrogen emission at moderately high redshifts. 
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1 Introduction 

In many astronomical imaging applications, images 
are taken at multiple wavelengths. Although the abil- 
ity to detect faint sources can be enhanced by stack- 
ing the images, a simple weighted linear combination 
produces a spectral bias dictated by the particular 
weighting function. We describe a detection algorithm 
which overcomes this limitation, and discuss its ap- 
plication to the Wide- field Infrared Survey Explorer 
(WISE) mission ijWright et aljl201(f ). 

Over a period of seven months from a precessing- 
polar orbit of the Earth, WISE surveyed the whole 
sky in four infrared bands with effective wavelengths 
of 3.4 ^m (Wl), 4.6 /im (W2), 12 /im (W3) and 22 
/im (W4) using an imaging array with a pixel size 
of 2.75", and spatial resolution (FWHM) of approx- 
imately 6" at the three shortest wavelengths and 12" 
at the longest. The WISE bands were chosen to opti- 
mize detection of cool brown dwarfs and luminous in- 
frared galaxies, but they are also well placed to study 
most objects in the universe. In general, the short 
wavelength bands are sensitive to starlight, while the 
long wavelength bands are sensitive to emission from 
the interstellar medium and fro m dust associa t ed wi th 
star formation (see, for example. I Jarrett et al.l l|201ll )'). 
We have performed source detection on the resulting 
stacks of four-band images using the Multiband DE- 
Tection (MDET) algorithm which is optimal for the 
detection of point sources in the presence of additive 
Gaussian noise. In the context of WISE, MDET rep- 
resents the initial detection step in source photometry. 
Its role is to produce a set of candidate sources which 
are then forwarded to a separate module for detailed 
parameter estimation consisting of source position, the 
flux at each band, the corresponding uncertainties, and 
various measures of the estimation quality. 

We discuss our experience with MDET using WISE 
data, and discuss its potential benefits for source de- 
tection with the Australian Square Kilometre Array 



Pathfinder ( ASKAP), whose design characteristics are 
discussed by[jphnston ct al. (2009|). In particular, we 
discuss its applicability to the search for neutral hy- 
drogen in various redshift ranges, as planned for key 
projectfQ WALLABY (Widefield ASKAP L-Band Legacy 
All-Sky Blind survey) and FLASH (The First Large 
Absorption Survey in HI). 

2 Theoretical Basis 
2.1 Measurement Model 

The starting point for the detection step is the mea- 
surement model for an isolated point sourceQ, assumed 
to be at location s and to have flux fx in the waveband 
denoted by index A; it can be expressed as: 

Px, = f\Hx(r\i - s) + b Xl + vxi (1) 

where p\i is the observed value of the ith pixel at sky 
location rxi, H\(y) is the point spread function (PSF) 
representing the response of a focal-plane pixel to a 
point source, b\i is the background sky level, and vx% is 
the noise, assumed to be a zero-mean Gaussian random 
process with covariance C v defined by: 

(C„)xi,X'i' = EvxiVx'i' — &XX'&ii'{Ov)xi + (Cb)\i,X'i' 

(2) 

where E is the expectation operator, (<T v )xi is the stan- 
dard deviation of measurement noise, assumed to be 
uncorrelated, and matrix Cb represents the covariance 
of the background. 

It is advantageous to estimate the background, bxi, 
ahead of time (using, for example, median filtering 
with a window size, W, appropriate to the characteris- 
tic spatial scale of background variations) and subtract 

1 http: / /www. atnf.csiro.au/SKA/ssps. html 
2 In Subsection 2.3 we will discuss the behavior in 
crowded fields, where this assumption is often violated 
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its contribution, so that the measurement model may 
be rewritten: 

PM = f\Hx(r\i - s) + v\i (3) 

If the observations are not sky background limited, 
or if the sky background is flat, the fact that a back- 
ground has been subtracted is not an issue. Otherwise, 
the analysis becomes complicated by the presence of 
spatial correlations in the residuals after subtraction. 
To minimize the effects of such correlations, W should 
be chosen comparable to the minimum spatial scale, 
L, of background fluctuations. Provided the correlated 
component of background residuals is not too large in 
comparison to the measurement noise, i.e., provided 
the diagnonal elements of C;, are smaller than (cr 2 )Ai, 
then the residuals of sky background subtraction can 
be lumped in with the measurement noise. Specifi- 
cally, the variance of subtraction residuals can then 
be treated as a "background confusion" term which is 
added in quadrature to (<T„) Ai . In the case of WISE, 
this was a good approximation since the residuals of 
background subtraction were dominated by low-level 
unresolved point sources and/or diffuse nebulosity on 
significantly larger scales than the sources of interest. 
The former component can be treated as uncorrelated 
noise for present purposes, while the latter component 
normally subtracts out cleanly. An example involving 
background subtraction in a region of heavy nebulos- 
ity is presented in Section 3. In subsequent discussion 
we therefore assume that the residuals of background 
subtraction can be treated in a similar fashion to un- 
correlated measurement noise, and denote the variance 
of combined noise as cr^i- 

There could conceivably be situations in which the 
above assumption would not be appropriate. For ex- 
ample, if background fluctuations were of such a scale 
as to be difficult to distinguish from genuine sources, 
then a more sophisticated approach would need to be 
taken. This would not invalidate the detection ap- 
proach to be described in the next section, however, 
because the treatment could be extended to the case of 
non-negligible background correlations in a straighfor- 
ward way. It would involve replacing the factor 1/afi, 
which appears in various summations, by {C^ )xi,X'i' i 
the appropriate summations would then be performed 
over the set of A, i, A', i' rather than simply A, i. 



2.2 Detection Algorithm 

Based on the measurement model expressed by Equa- 
tion (O, the source detection procedure involves com- 
paring the relative probabilities of the following two 
hypotheses at each location, s, within a predefined reg- 
ular grid of points on the sky: 

Hypothesis (A): s lies on blank sky at all wave- 
lengths 

Hypothesis (B): s represents the location of a source 
whose flux densities are the most probable values, de- 
noted by fx. 



2.2.1 Prior information on flux values: the 
positivity constraint 

To compare the above hypotheses requires knowledge 
of fx, which we obtain by maximizing the conditional 
probability, P(fjp), with respect to f, where f is a 
vector whose components are the set of fx, and p is a 
vector whose components are the set of pixel values, 
pxi, in the vicinity of s. The conditional probability 
itself is given by Bayes' rule, i.e., 

P(t\p) = P(p\f)P(f)/P{p) (4) 

where 

\nP(p\f) = -i V 4" [PXi ~ hH x (r Xi - s)] 2 (5) 
1 x,i 

and P(f ) represents our a priori knowledge about pos- 
sible flux values. We make no a priori assumptions 
about relative probabilities of spectral shapes of astro- 
physical objects, so our P(f ) is completely neutral on 
that point, i.e., we do not wish to introduce any color 
biases into the detector. However, one important piece 
of knowledge that we do have is that flux is positive. 
We can thus express P(f) as: 

P(f) = / COI1St - if h ^0 VA (Q\ 

* ' | otherwise ^ ' 

The remaining quantity, P(p), in ijl)). represents 
a normalization factor. The maximization of P(f\p) 
then yields: 




where the point spread function has been abbreviated 
to Hxi = Hx(rxi — s), and 8(x) represents a function 
which is equal to its argument if the latter is nonnega- 
tive and otherwise. The summations in (0) are over 
all pixels within a predefined neighborhood of s. 

2.2.2 Evaluating the relative probabilities 

With a further application of Bayes' rule, we can now 
express the probabilities of hypotheses (A) and (B), 
above, as: 

P(sky|p, m ) = P(sky, m ) ]^[ P(pA|sky, m )/P(px,m ) 

x 

(8) 

P(f\p, m) = P(f , m) J| P(p x \f, m)/P{px, m) (9) 

A 

where mo represents the sky-only model corresponding 
to hypothesis (A), and m represents the model corre- 
sponding to hypothesis (B), based on Equation (|3J. 

The likelihoods P(p|sky, mo) and P(p|/A,m) are 
given by: 

1 2 

lnP(pA|sky,m ) = --^^+const. (10) 

2 
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\nP(px\fx,m) 



\ E ~T^P X% ~ h H \if + const. 



(11) 

Assuming that we have no prior knowledge about 
the possible presence or absence of a source at s, all 
of the other factors in JH} and ([9]) may be regarded 
as constants for present purposes, and we can thus 
express the probability ratio (source/sky) as: 



In 



P{i\p,m) 
P(sky|p,m ) 



lJ2^J2^+ cowt - ( 12 ) 



Substituting for fx using Q, we can express this 
probability ratio as: 



In 



P(sky|p,m ) 
where <b(s) is defined as: 



i</>(s) 2 + const. 



E 



fl(E>A>L)gA(r Al -s) 
^(l/aDHxirxi-sy 



(13) 



(14) 



2.2.4 Summary of the detection procedure 

The multiwavelength detection algorithm as derived 
above consists of the following steps: 

1. Subtract a slowly varying background from the 
images at each of the individual wavelengths in 
order to "flatten" the sky. 

2. Calculate a spatial matched filter image at each 
individual wavelength. The result, obtained by 
cross-correlating the observed image with the 
PSF, optimizes the S/N of the point sources in 
the image at that wavelength. 

3. Divide each such matched filter image by a cor- 
responding uncertainty image representing the 
spatial variation in the standard deviation of lo- 
cal background noise. 

4. Set negative pixel values in the resulting image 
to zero. 

5. Form the quadrature sum of the clipped images. 

6. Threshold this image at the desired signal-to- 
noise level, Td, and find all local maxima in the 
thresholded image. 

A graphical illustration of the image combination 
procedure is given in Figure [T] 



in which we have replaced the values of the point spread 
function, Hxi, by their more explicit form. 

2.2.3 Criterion for source detection 

From (|13|1 . maxima in 0(s) correspond to maxima in 
the (source/sky) probability ratio, and hence an image 
formed by calculating <f)(s) over a regular grid of posi- 
tions, s, would be a suitable basis for optimal source 
detection. It is apparent from (|14[) that such an im- 
age represents a quadrature sum of matched filters at 
the individual wavelengths, with appropriate normal- 
ization. The noise properties of such an image can be 
assessed by expressing </>(s) in terms of the a posteriori 
variance of fx, given by: 




(S/N)a 



(S/N 



Hi 



(15) 



from which we obtain: 




(16) 



It is readily shown, from (|16p . that the standard 
deviation of <^>(s) is unity, i.e., <f)(s) itself is in units 
of standard deviations. Therefore, for a given detec- 
tion threshold Td [sigmas], the most likely locations 
of sources correspond to those for which (f)(s) > Td. 
The quantity Td represents our pre-defined detection 
threshold, which for the WISE Preliminary Release 
(jCutri et al.ll201lD was set at 7, i.e., all peaks above 
a signal to noise ratio of 7 were taken as candidate 
sources which were then passed to the photometry 
module for precise estimation of flux and position. 



Figure 1: Graphical illustration of the effect of 
combining matched filter images at multiple wave- 
lengths, three in this case. The axes in the fig- 
ure represent the pixel values corresponding to 
the source peak in the three individual bands, in 
units of the standard deviation of local measure- 
ment noise. When the corresponding matched fil- 
ter images are combined in quadrature, as per the 
MDET procedure, that pixel receives a value in- 
dicated by the length of the vector "Total." Since 
the signals in the individual bands are orthogonal, 
noise in individual bands has minimal effect on the 
resultant, and the latter is independent of which 
particular bands have the largest S/N, i.e., there 
is no bias towards any particular spectral shape. 



The MD ET detection algori thm is similar to one 
proposed by ISzalav et alj (1 19991 ) , often referred to as 
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the "chi squared" method. The central operation, in 
both cases, is a quadrature sum of matched filter im- 
ages. An important difference between the two pro- 
cedures, however, is the fact that in our technique 
we threshold the matched filter images at zero before 
squaring and combining, thus avoiding the contami- 
nating effect of squared negative values in the image 
sum. This is a direct consequence of our prior infor- 
mation concerning the positivity of intensity, imposed 
via our Bayesian framework using Equation It re- 
duces the background noise on the combined image by 
a factor of \/2, and therefore increases the sensitivity 
by the same factor. In applying the positivity con- 
straint, however, care must be exercised in the back- 
ground subtraction step, particularly in a confused re- 
gion in which the background may vary on relatively 
short spatial scales. 



2.3 Allowance for confusion 

In the above derivation, the assumption was made 
that the images of adjacent sources do not overlap. 
Although all matched filters are subject to this lim- 
itation, the effects of confusion can be can be more 
acute in the multiwavelength case. For example, in a 
WISE combined detection image, a star in Wl may 
become confused with an extended source such as a 
galaxy in W4. We overcome such effects by supple- 
menting the set of multiband detections with the re- 
sults of single-band detections, thereby maintaining 
the increased sensitivity of multiband detection while 
not missing any detections due to cross-band blend- 
ing effects. Most WISE sources, however, are detected 
in the multiband step; the number of supplementary 
single-band detections is typically only ~ 1% of the 
total, and many of those are simply the result of noise 
bumps. We include them for considerations of com- 
pleteness. 

The fact that the mathematical formalism on which 
MDET is based makes no explicit allowance for con- 
fusion is not a problem within the overall scheme of 
WISE source photometry, since MDET is simply the 
initial stage of a procedure in which the results are 
subsequently refined in a parameter estimation step. 
For example, if two closely spaced components are 
blended into a single peak in the MDET detection im- 
age, they are subsequently separated by the so-called 
"active d eblending" proced ure in profile-fitting pho- 
tometry ()Cutri et al.ll201lh . Such a situation is sig- 
naled by the presence of an elevated value of the re- 
duced chi squared, xl, of the maximum likelihood fit 
and indicates that an extra point-source component 
must be added to the model. Extended sources rep- 
resent another violation of the assumptions of the de- 
tection algorithm, although this was not a serious is- 
sue for WISE since most of the sources were spatially 
unresolved. Modification of the algorithm to opti- 
mize MDET for the detection of faint extended sources 
would require the use of a set of extended source tem- 
plates in Equation (|I4|I instead of the point spread 
function, H\{r). 



3 Examples of Application to 
WISE Data 

Figure shows three examples in which WISE images 
at four wavelength bands have been combined to pro- 
duce a detection image using the MDET procedure 
described above. These examples serve to illustrate 
several aspects of the procedure. 

The first example is of a young stellar object (a 
class I protostar candidate) in the L1689 starforming 
cloud of p Oph. It is visible in all four WISE bands 
against a background of diffuse emission. 

The second e xample is of an u l tra-co ol brown dwarf, 
using data from iMainzer et al ] l|201ll ). This object 
(WISEPC J04583. 90+643451.9) has an estimated tem- 
perature of 600 K and a spectral class T9. At temper- 
atures such as these, the spectrum is sharply peaked 
near 4.5 /j,m, corresponding to a relatively narrow "is- 
land" of low opacity between the heavy absorption due 
to such molecular components as water and methane. 
For this reason, the source is by far the most promi- 
nent in the W2 waveband of WISE. The Wl and W2 
filters were, in fact, optimized for the detection of this 
feature. 

The third example is of a Hyper-Luminous InfraRed 
Galaxy (HyLIRG), using data from Eisenhardt et al. 
(2011, in preparation). This is a very red object, and 
is thus brightest in W4. 

In all three cases, the images in the individual 
bands have been convolved with the respective PSFs 
to produce, in essence, optimal matched filter images 
at those wavelengths, so the corresponding S/N values 
are directly comparable with that of the combined (de- 
tection) image. The actual S/N values are presented 
in Table □ 

Table 1: Detection S/N for WISE observations of 
three selected objects 

Object Wl W2 W3 W4 Combined 

(3.4 (im) (4.6 ftm) (12 fim) (22 fim) 
YSO 23.67 65.07 33.27 25.89 101.89 
BD 9.33 67.48 2.38 0.03 68.16 

HyLIRG 2.11 1.12 22.69 24.35 33.15 



The first example (YSO) demonstrates the improve- 
ment in detectability that results from combining the 
images at all four wavelengths. As discussed above, 
one of the steps involved in this procedure is to sub- 
tract a slowly-varying sky background. The latter was 
estimated by median filtering using a moving window 
of size 21" x 21", chosen to be representative of the spa- 
tial scale of the background variations. Note that the 
combined S/N (101.9) exceeds the quadrature com- 
bination of the S/N at the individual bands (81.1); 
the additional improvement is due to the effect of the 
positivity constraint. 

Such a gain in S/N is not obtained in the second 
example (BD) since the source is detected primarily 
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Figure 2: Examples of the use of MDET for 
combining images at multiple wavelengths to in- 
crease detection sensitivity. Top row: Young stel- 
lar object (YSO); Middle row: Ultra-cool brown 
dwarf (BD); Bottom row: Hyper-Luminous In- 
fraRcd Galaxy (HyLIRG). Coadded images, each 
with 5' x 5' held of view, are shown for the four 
WISE wavebands, Wl (3.4 fim), W2 (4.6 fjm), W3 
(12 jitm), and W4 (22 /jm); these were combined 
to produce the optimal detection image (labeled 
"COMBINED") shown at the right of each row. 
For each of the three objects, all five images are 
presented on the same (linear) intensity scale in 
units of S/N, whose peak values are 102, 68, and 
33 for the YSO, BD, and LyLIRG, respectively. 
In each case, the object of interest is indicated by 
the blue arrow for the waveband with highest S/N 
(W2 for the YSO and brown dwarf; W4 for the 
HyLIRG). 
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in one band only. The example does, however, demon- 
strate that non-detection in the three "dropout" bands 
did not contaminate the detection. The third example 
(HyLIRG) represents a combination of both circum- 
stances, since the source is detected in only two of the 
four bands. Nevertheless, the two bands in which the 
source is detected have served to increase the source 
S/N in the combined image, while the two dropout 
bands have not appreciably contaminated the result. 

Overall, the use of MDET as the initial detection 
step has benefited WISE photometry in the following 
ways: 

1. The detection sensitivity has been increased for 
the reasons discussed above. This increase does 
not come at the expense of reliability, since the 
image combination process tends to average out 
the noise bumps. The inferred reliability meets 
the WISE Func tional Requireme nt value of 99.9% 
for SNR > 20 l|Cutri et al.ll201ll ). 

2. The fact that a single position is supplied for 
each source, even though the source is detected 
in multiple bands, facilitates simultaneous multi- 
band parameter estimation. The advantages of 
the latter are that no bandmerging is required 
(thus avoiding band-to-band matching ambigui- 
ties crowded fields) , and that dropout bands au- 
tomatically receive a flux upper limit. 

4 Application to ASKAP 

MDET can facilitate the optimal detection of faint 
sources by combining images along the frequency axis 
of an ASKAP data cube without introducing biases 
which favor one spectral shape over another. However, 
because of the limited frequency range of ASKAP (700- 
1800 MHz receiver range; 300 MHz correlator band- 
width), the benefits of the algorithm would be real- 
ized to a much larger extent for sources whose spec- 
tra contain narrow-band features than for sources with 
broadband spectra. As an example of the latter, con- 
sider a typical AGN whose flux density spectrum in 
the ASKAP frequency range is of the approximate 
form S v oc v^ 1 . Although the combining of images 
over all ASKAP channels will increase the S/N by ap- 
proximately the square root of the number of channels 
involved, the same will also be true for a uniformly- 
weighted linear combination of images of such a source. 
In fact, for a source with a similar spectral index, the 
improvement in S/N afforded by MDET over that 
obtained with the linear combination would be less 
than 4%. The situation is radically different for spec- 
troscopy, however, and we now discuss an important 
potential application. 

A major scientific goal for ASKAP and SKA is the 
detection of neutral hydrogen in distant galaxies, an 
important aspect of the study of galaxy formation an d 
evolution (| Johnston et al.ll2008l ; lRawlings et al.ll2004T ). 
The MDET procedure could be used to advantage in 
such searches by providing a spectrally neutral way 
to combine the images for all narrow-band channels 
within the frequency range corresponding to a given 
range of redshifts. As an example, we consider the 



combination of signals from a WALLABY data cube 
whose spectral sampling interval corresponds to 4 km 
s _1 . We suppose that somewhere in the redshift range 
z < 0.1 is the signal from a velocity-broadened HI 
line, a typical width for which might be ~ 400 km 
s _1 l|Koribalski wM). Since the achievable gain in de- 
tectability increases monotonically with the per-channel 
S/N of the spectral peaks of the source, it is desirable 
to carry out the detection step at a spectral resolution 
commensurate with the spectral structure of interest, 
so that some boxcar averaging of the spectral channels 
may be warranted. In the present example, combin- 
ing the raw channels in groups of 25 would result in 
an effective channel width of 100 km s _1 , so that our 
source signal would be spread over 4 such channels. 
For a peak S/N per channel ~ 10, an unweighted lin- 
ear combination of these channels would then produce 
a signal with S/N ~ 2, i.e., barely detectable, whereas 
the MDET procedure would result in S/N ~ 12. Of 
course, the peak signal from a properly-tuned spec- 
tral matched filteop would be even greater, but we 
would then be optimizing for a particular line width 
and shape, and be less sensitive to other potentially in- 
teresting structure. For example, if optimized for 100 
km s _1 structure, the spectral matched filter could do 
no better than S/N ~ 10. While this numerical ex- 
ample was illustrative only, we can make the general 
statement that by combining images from WALLABY 
data cubes using the MDET procedure, we obtain an 
optimal answer to the question: In which portions of 
the field of view are the observations inconsistent with 
random noise? 

With an appropriate sign change in Equation (|14|1 . 
the procedure could be applied to the detection of HI 
absorption, and therefore be of potential benefit to the 
processing of data from the FLASH project in the red- 
shift range 0.5 < z < 1.0. Ultimately, it is expected 
that the full SKA will enable detection out to z ~ 2.5 
(|Zwaanll200rJ l. At such high redshifts, the emission is 
too weak to detect the HI clouds of individual galaxies. 
The problem can be mitigated, however, by combining 
the HI images o f overlapping clouds of galaxies with 
known redshifts (|Kh andai et al. 20TlT) . This approach 
enabled lChang et al .1 lj2010T ) to detect neutral hydrogen 
out to z = 0.8. MDET has the potential for further 
improving this technique, since it provides a way of 
combining the images without prior knowledge of the 
redshifts of individual clouds. 

5 Conclusions 

MDET provides a procedure for combining the narrow- 
band images within a data cube in such a way as to 
increase optimally the detection signal to noise ratio 
without introducing a color bias. We have used it suc- 
cessfully in the WISE mission, and believe it will be of 
substantial benefit to ASKAP, particularly for sources 
with narrow-band spectral features. In particular, we 
suggest that it could aid in the detection of neutral 
hydrogen in distant galaxies. 

3 As distinct from the spatial matched filter which is an 
integral part of the MDET procedure 
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